﻿/*Copyright (C) 2010-2012 RAAFCOM B.V. All rights reserved.

This file is part of Velsto and may be used under the terms of:
The GNU General Public License version 3 or higher, as published by
the Free Software Foundation.

This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.*/
// ******************************************************************************************************
// <copyright file="Kriging.cs" company="Caroso.inc">
//     Copyright (c) Carlos Osorio All rights reserved.
// </copyright>
// ******************************************************************************************************
// The contents of this file are subject to the Mozilla Public License Version 1.1 (the "License"); 
// you may not use this file except in compliance with the License. You may obtain a copy of the License at 
// http:// www.mozilla.org/MPL/ 
// Software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTY OF 
// ANY KIND, either express or implied. See the License for the specificlanguage governing rights and 
// limitations under the License. 
// 
// The Original Code is this .mapWindow 
// 
// The Initial Developer of this version of the Original Code is Carlos Osorio carosocali@gmail.com
// 
// Contributor(s): (Open source contributors should list themselves and their modifications here). 
// Change Log: 
// Date           Changed By      Notes
// 22 Nov 2010    Carlos Osorio   Inital upload t
// 28 Apr 2012    Maurits van Ravestijn(RAAFCOM BV) Velsto version reading xyz file
// ******************************************************************************************************
namespace Velsto
{
    using System;
    using System.IO;
    using System.Collections;
    using MapWinGeoProc;

    /// <summary>
    /// Main class  for Kriging Interporlation
    /// </summary> 
    public class Kriging
    {
       #region propierties
        /// <summary> Number of points</summary> 
        public int npoints = 0;

        /// <summary>This is the Binding</summary> 
        private double[,] mlags;
        
        /// <summary>Matrix with the coordinates of all points</summary> 
        private double[,] cpoints;

        /// <summary> Distance max </summary> 
        private double maux = 0;

        /// <summary>Max value of Semivariance</summary> 
        private double mauxV = 0;

        /// <summary>Initial lag.</summary> 
        private double lag = 0;

        /// <summary>Model defintion</summary> 
        private IKrigingModel mmodel;
        
        /// <summary>Index points</summary>
        private QuadTree quadTree;

        /// <summary>Extention points xmin, ymin xmax, ymax</summary>
        private double[] extent = { double.MaxValue, double.MaxValue, double.MinValue, double.MinValue };

        /// <summary>
        ///  Min and Max values found in the interpolation
        /// </summary>
        //private double[] minmax = new double[] { double.MaxValue, double.MinValue };

        public  MapWinGIS.Grid outGrid;

        #endregion

        /// <summary>
        /// Initializes a new instance of the Kriging class
        /// </summary>
        /// <param name="shapeLayer">Shapefile type point</param>
        /// <param name="idField">Position of field that contains the Z value</param>
        public Kriging(MapWinGIS.Shapefile shapeLayer, int idField)
        {
            var modelLineal = new ModelLineal
                                  {
                                      C0 = 0, C1 = 1, Range = 1
                                  };
            this.CreatePoints(shapeLayer, idField);
        }

        public Kriging(MapWinGIS.Shapefile shpResid)
        {
            var modelLineal = new ModelLineal
            {
                C0 = 0,
                C1 = 1,
                Range = 1
            };
            this.CreatePoints(shpResid);
        }

        /// <summary>
        /// Initializes a new instance of the Kriging class
        /// </summary>
        /// <param name="shapeLayer">Shapefile (POINTS)</param>
        /// <param name="idField">Position of field that contains the Z value</param>
        /// <param name="model">This class provides the model that will be used</param>
        public Kriging(MapWinGIS.Shapefile shapeLayer, int idField, IKrigingModel model)
        {
            if (model != null)
            {
                this.mmodel = model;
            }

            this.CreatePoints(shapeLayer, idField);
            ////  GenerateMatrix(ref msemiVarTeo, ref inverMatrix, ref mdistances);
        }

        /// <summary>
        /// Initializes a new instance of the Kriging class
        /// </summary>
        /// <param name="resFile">File that contains xyz values</param>
        /// <param name="model">This class provides the model that will be used</param>
        public Kriging(MapWinGIS.Shapefile shpResid, IKrigingModel model)
        {
            if (model != null)
            {
                this.mmodel = model;
            }

            this.CreatePoints(shpResid);
            ////  GenerateMatrix(ref msemiVarTeo, ref inverMatrix, ref mdistances);
        }

        #region fields

        /// <summary>
        /// Gets or set the Min and Max values found in the interpolation
        /// </summary>
        //public double[] Minmax
        //{
        //    get { return minmax; }
        //    set { minmax = value; }
        //}

        /// <summary>Gets a value indicating number of points used by search</summary> 
        public int Np
        {
            get { return this.npoints; }
        }

        /// <summary>Gets a value indicating max Distance</summary> 
        public double MaxDistance
        {
            get { return this.maux; }
        }

        /// <summary>Gets a value indicating max variance</summary> 
        public double MaxVariance
        {
            get { return this.mauxV; }
        }

        /// <summary>Gets a value indicating number of point in the Shapefile</summary> 
        public int NumPoints
        {
            get { return this.npoints; }
        }

        /// <summary>Gets or sets the  length of the lag</summary> 
        public double ValueLag
        {
            get { return this.lag; }
            set { this.lag = value; }
        }

        #endregion

        /// <summary>
        /// This function creates the HDR of Gridfile
        /// </summary>
        /// <param name="inExtents"> Extension of grid</param>
        /// <param name="cellSize">Size cell of the grid</param>
        /// <param name="projection">Projection (the same that shapefile)</param>
        /// <param name="noDataValue">No value definition</param>
        /// <param name="outGridPath">Path of the output</param>
        /// <param name="outGrid">Name of the output grid</param>
        public static void CreateGridFromExtents(
            MapWinGIS.Extents inExtents, double cellSize, string projection, double noDataValue, string outGridPath, out MapWinGIS.Grid outGrid)
        {
            double height = Math.Abs(inExtents.yMax - inExtents.yMin);
            double width = Math.Abs(inExtents.xMax - inExtents.xMin);

            outGrid = new MapWinGIS.Grid();
            MapWinGIS.GridHeader hdr = new MapWinGIS.GridHeader();

            hdr.dX = cellSize;
            hdr.dY = cellSize;
            hdr.NodataValue = noDataValue;
            hdr.NumberRows = int.Parse(Math.Ceiling(height / cellSize).ToString());
            hdr.NumberCols = int.Parse(Math.Ceiling(width / cellSize).ToString());
            hdr.Projection = projection;
            hdr.XllCenter = inExtents.xMin + (0.5 * cellSize);
            hdr.YllCenter = inExtents.yMin + (0.5 * cellSize);
            outGrid.CreateNew(outGridPath, hdr, MapWinGIS.GridDataType.DoubleDataType, noDataValue, true, MapWinGIS.GridFileType.UseExtension, null);
        }

       /// <summary>
       /// This method captures the points of the shapefile
       /// </summary>
       /// <param name="shapeLayer">Shapefile (POINTS)</param>
       /// <param name="idField">Position of field that contains the Z value</param>
        public void CreatePoints(
            MapWinGIS.Shapefile shapeLayer, int idField)
        {
            this.npoints = shapeLayer.NumShapes;
            this.cpoints = new double[this.npoints, 3];
            this.quadTree = new QuadTree();

            for (int shp = 0; shp < shapeLayer.NumShapes; shp++)
            {
                MapWinGIS.Point pt = shapeLayer.QuickPoint(shp, 0);
                this.cpoints[shp, 0] = pt.x;
                this.cpoints[shp, 1] = pt.y;
                ////      this.cpoints[shp, 2] = (double)shapeLayer.get_CellValue(idField, shp);
                this.cpoints[shp, 2] = Convert.ToDouble(shapeLayer.get_CellValue(idField, shp));
                this.quadTree.Insert(pt.x, pt.y, shp);

                if (pt.y > this.extent[3])
                {
                    this.extent[3] = pt.y;
                }

                if (pt.x > this.extent[2])
                {
                    this.extent[2] = pt.x;
                }

                if (pt.x < this.extent[0])
                {
                    this.extent[0] = pt.x;
                }

                if (pt.y < this.extent[1])
                {
                    this.extent[1] = pt.y;
                }
            }

          ////  MessageBox.Show(extent[0].ToString() + " - " + extent[1].ToString() + "--" + extent[2].ToString() + " - " + extent[3].ToString());
            this.CalculateDistances();
        }

        // Utility Method to parse with both decimal separators
        public double parseDoubleString(string sd)
        {
            var currentCulture = System.Globalization.CultureInfo.InstalledUICulture;
            var numberFormat = (System.Globalization.NumberFormatInfo)
            currentCulture.NumberFormat.Clone();
            numberFormat.NumberDecimalSeparator = ".";
            double number;
            bool succeeded = double.TryParse(sd, System.Globalization.NumberStyles.Any, numberFormat, out number);
            if (succeeded)
            {
                return number;
            }
            else
            {
                numberFormat.NumberDecimalSeparator = ",";
                bool succeeded2 = double.TryParse(sd, System.Globalization.NumberStyles.Any, numberFormat, out number);
                if (succeeded2)
                {
                    return number;
                }
                else
                {
                    return -999.25;
                }
            }
        }






        /// <summary>
        /// This method captures the points of the Residual File
        /// </summary>
        /// <param name="resFile">File with residuals</param>
        public void CreatePoints(MapWinGIS.Shapefile shpResid)
        {
            int numpoints = 0;
            double x = 0.0, y = 0.0 , z = 0.0;
            for (int i3 = 0; i3 < shpResid.NumShapes; i3++)
            {
                z = (double)shpResid.get_CellValue(1, i3);
                if (z < Constants.SHP_NODATA)
                {
                    numpoints++;
                }
            }
            //this.cpoints = new double[shpResid.NumShapes, 3];
            //this.npoints = shpResid.NumShapes;
            this.cpoints = new double[numpoints, 3];
            this.npoints = numpoints;
            this.quadTree = new QuadTree();
            int i2 = 0;
            for (int i4 = 0; i4 < shpResid.NumShapes; i4++)
            {
                shpResid.get_Shape(i4).get_XY(0, ref x, ref y);
                z = (double)shpResid.get_CellValue(1, i4);
                if (z < Constants.SHP_NODATA)
                {
                    this.cpoints[i2, 0] = x;
                    this.cpoints[i2, 1] = y;
                    ////      this.cpoints[shp, 2] = (double)shapeLayer.get_CellValue(idField, shp);
                    this.cpoints[i2, 2] = z;
                    this.quadTree.Insert(x, y, i2);
                    i2++;
                    if (y > this.extent[3])
                    {
                        this.extent[3] = y;
                    }

                    if (x > this.extent[2])
                    {
                        this.extent[2] = x;
                    }

                    if (x < this.extent[0])
                    {
                        this.extent[0] = x;
                    }

                    if (y < this.extent[1])
                    {
                        this.extent[1] = y;
                    }
                }
            }

            ////  MessageBox.Show(extent[0].ToString() + " - " + extent[1].ToString() + "--" + extent[2].ToString() + " - " + extent[3].ToString());
            this.CalculateDistances();
        }

        /// <summary>
        /// Calculate the distastance between two points
        /// </summary>
        /// <param name="i">Position ith in the point array</param>
        /// <param name="j">Position jth in the point array</param>
        /// <returns>Return the distance</returns>
        public double Mdistances(
            int i, int j)
        {
          return Math.Sqrt(((this.cpoints[i, 0] - this.cpoints[j, 0]) * (this.cpoints[i, 0] - this.cpoints[j, 0])) + ((this.cpoints[i, 1] - this.cpoints[j, 1]) * (this.cpoints[i, 1] - this.cpoints[j, 1])));
        }

        /// <summary>
        /// Calculate Semivariances
        /// </summary>
        /// <param name="i">Position ith in the point array</param>
        /// <param name="j">Position jth in the point array</param>
        /// <returns>the semivariance</returns>
        public double MsemiVarEmp(
             int i, int j)
        {
            return 0.5 * ((this.cpoints[i, 2] - this.cpoints[j, 2]) * (this.cpoints[i, 2] - this.cpoints[j, 2]));
        }

        /// <summary>
        /// This Method generate the surface
        /// </summary>
        /// <param name="inExtents">Extents of the shapefile</param>
        /// <param name="cellSize">Resolution of surface</param>
        /// <param name="projection">Projection shapefile</param>
        /// <param name="noDataValue">Value Null</param>
        /// <param name="outGridPath">Name of the file that will be created</param>
        /// <param name="nPn">Number of points that will be used to interpolate each cell </param>
        public void Surface(
            MapWinGIS.Extents inExtents, double cellSize, string projection, double noDataValue, string outGridPath, int nPn)
        {
           // this.CalculateDistances();
           // this.GenerateMatrix(ref this.mSemiVarTeo, ref this.inverMatrix, ref this.mDistances);

            DataManagement.DeleteGrid(ref outGridPath);
            ////  NoDataValue = -32768;
            CreateGridFromExtents(inExtents, cellSize, projection, noDataValue, outGridPath, out outGrid);

            int nr = outGrid.Header.NumberRows;
            int nc = outGrid.Header.NumberCols;

            double projX, projY;
            double[] cellPoint;
            for (int row = 0; row < nr; row++)
            {
                for (int col = 0; col < nc; col++)
                {
                    outGrid.CellToProj(col, row, out projX, out projY);
                    cellPoint = new double[2];
                    cellPoint[0] = projX;
                    cellPoint[1] = projY;
                    outGrid.set_Value(col, row, this.Interpolate(nPn, projX, projY));
                }
            }

            outGrid.Save(outGridPath, MapWinGIS.GridFileType.UseExtension, null);
           // outGrid.Close();
        }

        /// <summary>
        /// CalculateLags
        /// This method calculates the semivariance taking into account the lag
        /// </summary>
        /// <param name="lag">Distance. It depends on the average distances amogns the points</param>
        /// <param name="nlag">Number of Lags. Typically is 12.</param>
        /// <param name="gdist">Get distaces in an array</param>
        /// <param name="gsemi">Get semivariances in an array </param>
        public void CalculateLags(
            double lag, int nlag, out ArrayList gdist, out ArrayList gsemi)
        {
            ////0.min 1.Max 2. Dist sum, 3. Semi sum, 4 number of points
            this.mlags = new double[nlag, 5]; 
            gdist = new ArrayList();
            gsemi = new ArrayList();
            double auxR = 0;
            for (int i = 0; i < nlag; i++)
            {
                this.mlags[i, 0] = auxR;
                this.mlags[i, 1] = auxR + lag;
                this.mlags[i, 2] = 0;
                this.mlags[i, 3] = 0;
                this.mlags[i, 4] = 0;
                auxR += lag;
            }

            for (int i = 0; i < this.npoints; i++)
            {
                for (int j = i; j < this.npoints; j++)
                {
                    if (i != j)
                    {
                        double tdist = 0;
                        double tsemi = 0;
                        tdist = this.Mdistances(i, j);
                        tsemi = this.MsemiVarEmp(i, j);
                        for (int k = 0; k < nlag; k++)
                        {
                            if ((tdist >= this.mlags[k, 0]) && (tdist < this.mlags[k, 1]))
                            {
                                this.mlags[k, 2] += tdist;
                                this.mlags[k, 3] += tsemi;
                                this.mlags[k, 4] += 1;
                            }
                        }
                    }
                }
            }

            for (int i = 0; i < nlag; i++)
            {
                if (this.mlags[i, 4] > 0)
                {
                    this.mlags[i, 2] = this.mlags[i, 2] / this.mlags[i, 4];
                    this.mlags[i, 3] = this.mlags[i, 3] / this.mlags[i, 4];
                    gdist.Add(this.mlags[i, 2]);
                    gsemi.Add(this.mlags[i, 3]);
                }
                else
                {
                    this.mlags[i, 2] = this.mlags[i, 0] + (0.5 * lag);
                }
            }
        }

        /// <summary>
        /// Generate the arrays of distances and semivariances for drawing the theoretical semivariance
        /// </summary>
        /// <param name="gdistF">Array distance</param>
        /// <param name="gsemiT">Array semivariances</param>
        public void CalculateModelGraphic(
            out ArrayList gdistF, out ArrayList gsemiT)
        {
            double div = this.maux / 30;
            gdistF = new ArrayList();
            gsemiT = new ArrayList();
            for (int i = 0; i < 30; i++)
            {
                gdistF.Add(div * i);
                gsemiT.Add(this.mmodel.GetValue(div * i));
            }
        }
        
        /// <summary>
        /// This method generates an ArraList with distances for drawing the semivariogram
        /// </summary>
        /// <returns> An array with the distances</returns>
        public ArrayList Distances()
        {
            ArrayList mdist = new ArrayList();
            for (int i = 0; i < this.npoints; i++)
            {
                for (int j = i; j < this.npoints; j++)
                {
                    if (j != i)
                    {
                        mdist.Add(this.Mdistances(i, j));
                    }
                }
            }

            return mdist;
        }
        
        /// <summary>
        /// This method generates an ArraList with semivariances for drawing the semivariogram
        /// </summary>
        /// <returns>An array with all semivariances</returns>
        public ArrayList SemiVariances()
        {
            ArrayList msemi = new ArrayList();
            for (int i = 0; i < this.npoints; i++)
            {
                for (int j = i; j < this.npoints; j++)
                {
                    if (j != i)
                    {
                        msemi.Add(this.MsemiVarEmp(i, j));
                    }
                }
            }

            return msemi;
        }

        /// <summary>
        /// This method gets the value depending on the number of points to search nPn, and the position of the point
        /// </summary>
        /// <param name="nPn">Number of points to execute the interpolation</param>
        /// <param name="x">Coordinate X of the point to calculate</param>
        /// <param name="y">Coordinate Y of the point to calculate</param>
        /// <returns>Value found after interpolation</returns>
        public double Interpolate(
            int nPn, double x, double y)
        {
            double value = 0;
            double[] vectorG = new double[nPn + 1];
            double[] weights = new double[nPn + 1];
            vectorG[nPn] = 1.0;

            double[,] cpointsF;
            double[,] mdistF = new double[nPn + 1, nPn + 1];
            double[,] msemEF = new double[nPn + 1, nPn + 1];
            double[,] msemTF = new double[nPn + 1, nPn + 1];
            double[,] inverMatrixF = new double[nPn + 1, nPn + 1];

            this.ListPointAvailables(out cpointsF, nPn, x, y);
            this.CalculateDistancesAround(nPn, ref cpointsF, ref mdistF, ref msemEF);
            this.GenerateMatrixProximity(nPn, ref msemTF, ref inverMatrixF, ref mdistF);

            if (inverMatrixF == null)
                return double.NaN;


            for (int i = 0; i < nPn; i++)
            {
                vectorG[i] = this.mmodel.GetValue(Math.Sqrt(((cpointsF[i, 0] - x) * (cpointsF[i, 0] - x)) +
                        ((cpointsF[i, 1] - y) * (cpointsF[i, 1] - y))));
            }

            for (int i = 0; i < nPn + 1; i++)
            {
                weights[i] = 0;
                for (int j = 0; j < nPn + 1; j++)
                {
                    weights[i] += inverMatrixF[i, j] * vectorG[j];
                }
            }

            for (int i = 0; i < nPn; i++)
            {
                value += weights[i] * cpointsF[i, 2];
            }



            //if (this.minmax[0] > value)
            //{
            //    this.minmax[0] = value;
            //}

            //if (this.minmax[1] < value)
            //{
            //    this.minmax[1] = value;
            //}


            return value;
        }
 
        /// <summary>
        /// Searching the points around of a position
        /// </summary>
        /// <param name="cpointsF">Point founds</param>
        /// <param name="nPn">Number of point searched</param>
        /// <param name="x">Coordinate X of the point to calculate</param>
        /// <param name="y">Coordinate Y of the point to calculate</param>
        public void ListPointAvailables(
            out double[,] cpointsF, int nPn, double x, double y)
        {
            double mult = 1d;
            double searchDistance = (this.extent[2] - this.extent[0]) * 0.05;

            do
            {   
                this.quadTree.Clean();
                double xmin = x - (searchDistance * mult);
                double ymin = y - (searchDistance * mult);
                double xmax = x + (searchDistance * mult);
                double ymax = y + (searchDistance * mult);
                Interval intX = new Interval(xmin, xmax);
                Interval intY = new Interval(ymin, ymax);
                Interval2D rect = new Interval2D(intX, intY);
                this.quadTree.Query2D(rect);
                mult++;
            }
            while (this.quadTree.GetData().Count <= nPn);

            double[,] cpoint = new double[this.quadTree.GetData().Count, 3];
       
             int ic = 0;
             int[] idPoint = new int[this.quadTree.GetData().Count];
             foreach (object obj in this.quadTree.GetData())
             {
                 cpoint[ic, 0] = this.cpoints[(int)obj, 0];
                 cpoint[ic, 1] = this.cpoints[(int)obj, 1];
                 cpoint[ic, 2] = this.cpoints[(int)obj, 2];
                 idPoint[ic] = (int)obj;
                 ic++;
             }

            cpointsF = new double[nPn, 3];
            double[] mdist = new double[ic];

            for (int i = 0; i < ic; i++)
            {
                mdist[i] = Math.Sqrt(((cpoint[i, 0] - x) * (cpoint[i, 0] - x)) +
                        ((cpoint[i, 1] - y) * (cpoint[i, 1] - y)));
            }

            this.OrderPoints(ref mdist, ref idPoint, ic);

            for (int i = 0; i < nPn; i++)
            {
                cpointsF[i, 0] = this.cpoints[idPoint[i], 0];
                cpointsF[i, 1] = this.cpoints[idPoint[i], 1];
                cpointsF[i, 2] = this.cpoints[idPoint[i], 2];
            }
        }

        /// <summary>
        /// This method calculates the distance between points that are close to a point
        /// </summary>
        /// <param name="npoints">Number of points around</param>
        /// <param name="cpoints">Matrix with coordinates</param>
        /// <param name="mdistances">Distances between points</param>
        /// <param name="msemiVarEmp">Semivariance of point selected</param>
        private void CalculateDistancesAround(
            int npoints, ref double[,] cpoints, ref double[,] mdistances, ref double[,] msemiVarEmp)
        {
            for (int i = 0; i < npoints; i++)
            {
                for (int j = 0; j < npoints; j++)
                {
                    mdistances[i, j] = Math.Sqrt(((cpoints[i, 0] - cpoints[j, 0]) * (cpoints[i, 0] - cpoints[j, 0])) +
                         ((cpoints[i, 1] - cpoints[j, 1]) * (cpoints[i, 1] - cpoints[j, 1])));
                    msemiVarEmp[i, j] = 0.5 * ((cpoints[i, 2] - cpoints[j, 2]) * (cpoints[i, 2] - cpoints[j, 2]));

                    if (mdistances[i, j] > this.maux)
                    {
                        this.maux = mdistances[i, j];
                    }

                    if (msemiVarEmp[i, j] > this.mauxV)
                    {
                        this.mauxV = msemiVarEmp[i, j];
                    }
                }
            }
        }

        /// <summary>
        /// This method calculate the distances between all points
        /// </summary>
        private void CalculateDistances()
        {
            double mdis = 0d;
            double msem = 0d;
            for (int i = 0; i < this.npoints; i++)
            {
                for (int j = 0; j < this.npoints; j++)
                {
                    mdis = this.Mdistances(i, j);
                    msem = this.MsemiVarEmp(i, j);

                    if (mdis > this.maux)
                    {
                        this.maux = mdis;
                    }

                    if (msem > this.mauxV)
                    {
                        this.mauxV = msem;
                    }
                }
            }
        }

        /// <summary>
        /// Generate initial matrix for interpolation using the point around of unknown value
        /// </summary>
        /// <param name="nPn">Number of point used</param>
        /// <param name="msemiVarTeo">Matrix semivariances</param>
        /// <param name="inva">Matrix inverse</param>
        /// <param name="mdistances">Matrix of distaces</param>
        private void GenerateMatrixProximity(
            int nPn, ref double[,] msemiVarTeo, ref double[,] inva, ref double[,] mdistances)
        {
            for (int i = 0; i <= nPn; i++)
            {
                msemiVarTeo[i, nPn] = 1.0;
                inva[i, nPn] = 1.0;
            }

            for (int i = 0; i <= nPn; i++)
            {
                msemiVarTeo[nPn, i] = 1.0;
                inva[nPn, i] = 1.0;
            }

            msemiVarTeo[nPn, nPn] = 0;
            inva[nPn, nPn] = 0;
            for (int i = 0; i < nPn; i++)
            {
                for (int j = 0; j < nPn; j++)
                {
                    msemiVarTeo[i, j] = this.mmodel.GetValue(mdistances[i, j]);
                    inva[i, j] = msemiVarTeo[i, j];
                }
            }

            Matriz matriz = new Matriz(inva);
            Matriz inva1=null;
            try{
            inva1 = matriz.Inverse();
            }catch
            {
             inva= null;
             return;
            }
            
            inva= inva1.ToArray();

           ////int info = 0;
           //// matinv.matinvreport rep = new matinv.matinvreport();
           //// matinv.rmatrixinverse(ref inva, nPn + 1, ref info, ref rep);
        }
       
        /// <summary>
        /// Sort a list of points depending on distances to unknown value 
        /// </summary>
        /// <param name="a">Distance vector</param>
        /// <param name="b">Id point vector</param>
        /// <param name="x">Number of points</param>
        private void OrderPoints(
            ref double[] a, ref int[] b, int x)
        {
            int i;
            int j;
            double temp;
            int tempb;
            for (i = x - 1; i >= 0; i--)
            {
                for (j = 1; j <= i; j++)
                {
                    if (a[j - 1] > a[j])
                    {
                        temp = a[j - 1];
                        tempb = b[j - 1];
                        a[j - 1] = a[j];
                        a[j] = temp;

                        b[j - 1] = b[j];
                        b[j] = tempb;
                    }
                }
            }
        }
     }
}
